ggvote_plot_mfx3 <- function(ames_data,
                             desc_data,
                             model_ids,
                             add_descriptives = TRUE,
                             d = "cmr_arm_cwu",
                             x1 = "Homeownership",
                             x2 = "income",
                             x3 = "region",
                             x2_continuous = TRUE,
                             subset_x1 = NULL,
                             ylim = rep(NA_real_, 2L),
                             ytitle_left = "Marginal Effect",
                             ytitle_right = "",
                             xtitle = "",
                             font_size = 15L) {
  
  ## Data
  if (x2_continuous) {
    data <- ames_data %>% 
      dplyr::filter(model_id == model_ids[1]) %>%
      dplyr::bind_rows(
        ames_data %>%
          dplyr::filter(model_id == model_ids[2])
      ) %>%
      dplyr::filter(factor == d) %>%
      dplyr::mutate(Homeownership = ifelse(model_id == model_ids[1],
                                           "Renters",
                                           "Owners"))
    # x2_breaks_1 <- sort(unique(data[data$model_id == model_ids[1], x2]))
    # x2_breaks_2 <- sort(unique(data[data$model_id == model_ids[2], x2]))
    # x2_breaks_1_mod <- c(
    #   x2_breaks_1[1] - abs(diff(x2_breaks_1[1:2])),
    #   x2_breaks_1,
    #   x2_breaks_1[length(x2_breaks_1)] + 
    #     abs(diff(x2_breaks_1[(length(x2_breaks_1) - 1):length(x2_breaks_1)]))
    # )
    # x2_breaks_2_mod <- c(
    #   x2_breaks_2[1] - abs(diff(x2_breaks_2[1:2])),
    #   x2_breaks_2, 
    #   x2_breaks_2[length(x2_breaks_2)] + 
    #     abs(diff(x2_breaks_2[(length(x2_breaks_2) - 1):length(x2_breaks_2)]))
    # )
    desc <- desc_data %>%
      dplyr::filter(model_id %in% model_ids) %>%
      dplyr::filter(variable == x2) %>%
      group_by(model_id, !!as.name(x3)) %>%
      dplyr::mutate(
        x2_lower = ifelse(row_number() < n(),
                          breaks,
                          NA_real_
        ),
        x2_upper = ifelse(row_number() < n(),
                          dplyr::lead(breaks),
                          NA_real_
        )
      ) %>%
      dplyr::filter(not(is.na(x2_upper))) %>%
      dplyr::mutate(Homeownership = ifelse(model_id == model_ids[1],
                                           "Renters",
                                           "Owners"))
  } else {
    data <- ames_data %>% 
      dplyr::filter(model_id == model_ids[1]) %>%
      dplyr::bind_rows(
        ames_data %>%
          dplyr::filter(model_id == model_ids[2])
      ) %>%
      dplyr::filter(factor == d) %>%
      dplyr::left_join(
        desc_data %>%
          dplyr::filter(model_id %in% model_ids) %>%
          dplyr::filter(variable == x2) %>%
          dplyr::rename(!!as.name(x2) := levels) %>%
          dplyr::select(model_id, counts, all_of(x2), all_of(x3)),
        by = c("model_id", x2, x3)
      ) %>%
      dplyr::mutate(Homeownership = ifelse(model_id == model_ids[1],
                                           "Renters",
                                           "Owners")) %>% 
      dplyr::ungroup() 
  }
  
  ## Subset by x1
  if (not(is.null(subset_x1))) {
    if (all(is.na(ylim))) {
      ylim <- c(min(data$lower, na.rm = TRUE),
                max(data$upper, na.rm = TRUE))
    }
    
    data <- data %>%
      filter(!!as.name(x1) == subset_x1)
    
    if (x2_continuous) {
      desc <- desc %>%
        filter(!!as.name(x1) == subset_x1)
    }
  }
  
  
  ## Marginal effect
  vote_plot <- data %>%
    ggplot(aes(x = !!as.name(x2),
               y = AME)) +
    geom_point(position = position_dodge(width = 1)) +
    geom_errorbar(aes(ymin = lower,
                      ymax = upper),
                  position = position_dodge(width = 1),
                  width = 0) +
    geom_hline(yintercept = 0, col = "red") +
    theme(panel.spacing = unit(0.75, "lines"))  +
    theme(text = element_text(size = font_size)) +
    ylim(ylim) +
    xlab(xtitle) +
    scale_y_continuous(
      name = ytitle_left,
      sec.axis = sec_axis(trans = ~., 
                          name = ytitle_right,
                          labels = NULL,
                          breaks = NULL)
    )
  
  if (is.null(subset_x1)) {
    vote_plot <- vote_plot +
      facet_grid(. ~ get(eval(x3)) ~ get(eval(x1)))
  } else {
    vote_plot <- vote_plot +
      facet_grid(. ~ get(eval(x1)) ~ get(eval(x3)))
  }
  
  ## Descriptives
  if (add_descriptives) {
    if (is.null(subset_x1)) {
      effective_ymin <-
        ggplot_build(vote_plot)$layout$panel_scales_y[[1]]$range$range[1]
      effective_ymax <-
        ggplot_build(vote_plot)$layout$panel_scales_y[[1]]$range$range[2]
    } else {
      effective_ymin <- ylim[1]
      effective_ymax <- ylim[2]
    }
    
    if (x2_continuous) {
      vote_plot <- vote_plot +
        geom_rect(
          data = desc %>%
            dplyr::rowwise() %>%
            dplyr::mutate(x2_mid = mean(c(
              x2_lower, x2_upper
            ))) %>%
            dplyr::ungroup() %>%
            dplyr::mutate(
              counts = dplyr::coalesce(counts, 0L),!!as.name(x2) := NA_real_,
              AME = NA_real_
            ),
          aes(
            xmin = x2_lower,
            xmax = x2_upper,
            ymin = effective_ymin,
            ymax = ((counts / max(counts)) *
                      abs(effective_ymax - effective_ymin)) +
              effective_ymin
          ),
          alpha = .15,
          color = adjustcolor("gray20", alpha.f = .25)
        ) +
        xlim(range(data[[x2]]))
    } else {
      vote_plot <- vote_plot +
        geom_segment(
          aes(
            x = !!as.name(x2),
            xend = !!as.name(x2),
            y = effective_ymin,
            yend = ((counts / max(counts)) *
                      abs(effective_ymax - effective_ymin)) +
              effective_ymin
          ),
          size = 25,
          alpha = .1
        )
    }
  }
  
  ## Return
  return(vote_plot)
}
